arXiv:1506.05708vl [math.NA] 18Jun2015 


A weak Local Linearization scheme for stochastic differential 
equations with multiplicative noise 

J.C. Jimenez * C. Moral M. Selva^ 

June 19, 2015 


Abstract 

In this paper, a weak Local Linearization scheme for Stochastic Differential Equations (SDEs) with 
multiplicative noise is introduced. First, for a time discretization, the solution of the SDE is locally 
approximated by the solution of the piecewise linear SDE that results from the Local Linearization 
strategy. The weak numerical scheme is then defined as a sequence of random vectors whose first 
moments coincide with those of the piecewise linear SDE on the time discretization. The rate of 
convergence is derived and numerical simulations are presented for illustrating the performance of the 
scheme. 


1 Introduction 

During 30 years the class of local linearization integrators has been developed for different types of 
deterministic and random differential equations. The essential principle of such integration methods 
is the piecewise linearization of the given differential equation to obtain consecutive linear equations 
that are explicitly solved at each time step. This general approach has worked well for the classes of 
ordinary, delay, random and stochastic differential equations with additive noise. Key element of such 
success is the use of explicit solutions or suitable approximations for the resulting linear differential 
equations. Precisely, the absence of explicit solution or adequate approximation for linear Stochastic 
Differential Equations (SDEs) with multiplicative noise is the main reason of the limited application of 
the Local Linearization approach to nonlinear SDEs with multiplicative noise. For these equations, the 
available local linearization integrators are of two types: the introduced in [2] for scalar equations and 
the considered in [nmnn]. The former uses the explicit solution of the scalar linear equations with 
multiplicative noise, while the latter employs the solution of the linear equation with additive noise that 
locally approximates the nonlinear equation. 

Directly related to the development of the local linearization integrators is the concept of Local Linear 
approximations (see, e.g., [SJ[71[n]). These approximations to the solution of the differential equations 
are defined as the continuous time solution of the piecewise linear equations associated to the Local 
Linearization method. These continuous approximations have played a fundamental role for studying the 
convergence, stability and dynamics of the local linearization integrators for all the classes of differential 
equations mentioned above with the exception of the SDEs with multiplicative noise. For this last class of 
equations, the Local Linear approximations have only been used for constructing piecewise approximations 
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to the mean and variance of the states in the framework of continuous-discrete filtering problems (see 

El)- 

The purpose of this work is to construct a weak Local Linearization integrator for SDEs with mul¬ 
tiplicative noise based on suitable weak approximation to the solution of piecewise linear SDEs with 
multiplicative noise. For this, we cross two ideas: 1) as in [2, the use of the Local Linear approximations 
for constructing piecewise approximations to the mean and variance of the SDEs with multiplicative noise; 
and 2 ) as in [3], at each integration step, the generation of a random vector with the mean and variance 
of the Local Linear approximation at this integration time. For implementing this, new formulas recently 
obtained in for the mean and variance of the solution of linear SDEs with multiplicative noise are 
used, which are computationally more efficient than those formerly proposed in BM- Notice that this 
integration approach is conceptually different to that usually employed for designing weak integrators for 
SDEs. Typically, these integrators are derived from a truncated Ito-Taylor expansion of the equation’s 
solution at each integration step, and include the generation of random variables with moments equal to 
those of the involved multiple Ito integrals uniiii]. 

The paper is organized as follows. After some basic notations in Section 2, the new Local Linearization 
integrator is introduced in Section 3. Its rate of convergence is derived in Section 4 and, in the last section, 
numerical simulations are presented in order to illustrate the performance of the numerical integrator. 


2 Basic notations 

Let us consider the SDE with multiplicative noise 

pt fn pt 

Xt = Xt,+ f{s,Xs)ds + Y, 9^{s,Xs)dWl Vte[to,T], (1) 

^ to k—1 

where /, : [toi T] x K'’* —>■ K'’* are smooth functions, W ^,..., IT’” are independent Wiener processes on a 

filtered complete probability space ^11, S’, (T?i)t>to ’ adapted K'^-valued stochastic process. 

In addition, let us assume the usual conditions for the existence and uniqueness of a weak solution of Q 
with bounded moments (see, e.g., [10]). 

Throughout this paper, we consider the time discretization to = tq < ti < ■ ■ ■ < tn = T with 
T„+i — Tn < A for all n = 0,..., iV — 1 and A > 0. We use the same symbol K (•) (resp., K) for 
different positive increasing functions (resp., positive real numbers) having the common property to be 
independent of (t/j)^,^^ Moreover, stands for the transpose of the matrix A, and I-] denotes the 
Euclidean norm for vectors or the Frobenious norm for matrices. By Cp we mean the collection 

of all £-times continuously differentiable functions 5 : —)■ M such that g and all its partial derivatives 
of orders 1,2,... ,i have at most polynomial growth. 


3 Numerical method 

Suppose that Zn ~ Xr„ with n = 0,..., N — 1. 
yields 

/ (t,x) « (t„,z„) -h — 

whenever x ~ Zn and t Ki Tn. Therefore 


Set g^ = f ■ Taking the first-order Taylor expansion of g^ 
dg^ 

{Tn, Zn) {X - Zn) + (t„, Z„) {t - Tn) 

at 


Xt 



{B'^Xs + (s)) dW^ 


Vt € \Tn, , 


2 


( 2 ) 


with = s, B^ = (r„, z„) and 


dg'^ 


dg^ 


(^) — 9 ^n'} {Ttij ^n) ('^n^ ^n) ^n) ■ 

This follows that, for all t G [Tyi,r^_^i], Xt can be approximated by 

m „t 

Yt = Zr, + Y, [B^Ys + (s)) dWl yt G [r„, t„+i] 

k=0 


(3) 


which is the first order Local Linear approximation of Xt used in m- Hence, E(/) (hir„+i) 

for any smooth function </), and so ^t„+i might be weakly approximated by a random variable Zn+i such 
that the first moments of Zn+i — Zn be similar to those of Lir„+i ~ Zn- This leads us to the following Local 
Linearization scheme. 

Scheme 1 . Let , 77 ™, • • • , • ■ •, ??]v-i i-i-d. symmetric random variables having variance 1 

and finite moments of any order. For a given Zq, we define recursively jy by 


Zn+l — Ln Y ^^n ('^n+l) /^n (’^n+l) Lm ('Tn-i-i) T^ti, 

where gn = (Vnj ■ ■ • 1 dn (t), an {t) satisfy the linear differential equations 


(4) 


Ln{t) = Zn + J (B° (s) + 6 ° (s)) ds Vf G [r„,T„+i], 

(^) — ZnZn T Ln (s, (^)) ds yt G [tVi, 771 - 1 - 1 ] ■ 


Here 


(5) 

( 6 ) 


£„ (s, a) = a (H°) ^ +B°a'^ + /x„ (s) (&° (s))^ + &° (s) 77]!; (s) 

m 

+ ^ gn (s) (s))^ + hi (s) 77T (s) (sfc)T + (5) (,))T j ^ 

A.-1 

Remark 3.1. From ^ it follows that yin ijn+i) is the expected valued ofYr givenYr = Zn- Moreover, 
implies 

an {Tn+l) = E = Zn) ■ 

Remark 3.2. By construction, Scheme^^preserves the mean-square stability property that the solution of 
the linear equation dXt = XtFh^'^tYh^''^)dWf might have. For instance, if the trivial solution 

of the homogenous equation dXt = B^XtdWf is mean-square asymptotically stable. Scheme 

inherits this property. 

Remark 3.3. A key point in the implementation of Scheme [7] is the evaluation of just one matrix 
exponential for computing pin (t'k+i) find an (tVi-i-i) at each time step. Indeed, from Theorem 2 in 

Pn (Tn-ei) = Zn+ ( 7 ) 
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and 


( 8 ) 


vec{an (r„+i)) = 

where the matrices M.n, ^i, ^2 o^'^d the vector Un are given by 


■ A 

B 5 

Bi 

B 3 

B2 

Bi ' 


vec{znzl) 

0 

c 

Xd+2 

0 

0 

0 


0 

0 

0 

c 

0 

0 

0 


r 

0 

0 

0 

0 

2 

0 

, Ufi — 

0 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

0 

0 

0 
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^2 = [ 0dx{d^+d+2) Id 0dx5 ] and Ci = [ Xj2 Od'^x{2d+7) ], with matrices A, Bi, C and r defined by 


m 


\ B° 

n 

bO,l 

n 

B°zu + by' 


r\ 


c = 

0 

0 

1 

, r = 

^(d+l)xl 

1 

k^l 


0 

0 

0 




Bi = wec(^i) + jd^Zn, B2 = vec{P2) + PbZn, B3 = vecifi^), Bi = PaC, and B^ = Here C = 

[ Xd 0dx2 ] and 


k^l k^l k^l 

m m 

/?4 = © 6^ + E © &n’°, ^5 = 6°’' © + E ' © ^^ + S" © bt'\ 


fc=l 


/c=l 


being and b^^ defined via Q) as b^’^+b^’^ {s — t„) = b^ (s). The symbols vec, 0 and 0 denote the vec- 
torization operator, the Kronecker sum and the Kronecker product, respectively. Xd is the d— dimensional 
identity matrix. The matrix exponential in 0 and 0 can be efficiently computed via the Fade method 
with scaling and squaring strategy or via the Krylov subspace method in the case of large system of SDKs 
(see, e.g., m)- For autonomous equations or for equations with additive noise, the exponential matrix 
in 0 and & reduces to simpler forms m- 

Remark 3.4. For SDEs with additive noise. Scheme^ reduces to the weak order -1 Local Linearization 
scheme introduced in m- 


4 Rate of convergence 

Next theorem establishes the linear rate of weak convergence of Scheme when the drift and diffusion 
coefficients are smooth enough. 

Hypothesis 1. For any k = 0 ,... ,m we have g^ G Cp ([to, T] x . Moreover, 

\9^ {bix)\ < K {l + \x\) and 


dg^ 

Idt 


{t,x) 


dg^ 

dx 


{t,x) 


< K 


( 9 ) 


for all t G [toj T] and x gM.'^. 
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Theorem 4 . 1 . In addition to Hypothesis^ suppose that Xt^ has finite moments of any order and that 
for all (j) G Cp 

\Efi{Xt„)-Efi{zo)\<KA. 


'4 tmd 


Then, for all G Cp 


where Zm is given hy Scheme\^ 


\Efi{XT)-Ecj){zN)\<K{T)A, 


Theorem 4.1 is a straightforward result of Theorem 14 . 5.2 in [TO] and the two following Lemmata. 


Lemma 4 . 1 . Under the assumptions of Theorem \ 4 - 1 \ for any q > 1 we have 

\2q 


E ( J < iL (T) (l + E (kol"')) 


and 


E (\Zn+l - <K{T) (Tfc+i - TkY (l + 

for all n = 0 ,..., N — 1 . 

Proof. From Hypothesis it follows that | | < K and 

|6^(s)| <if(r)(i + |z„|) 


( 10 ) 

( 11 ) 


( 12 ) 


for all n = 0,..., — 1, k = 0,... ,m and s G [t„, r„+i]. Then, combining Gronwall’s lemma with ([^ 

gives 

\p-n is)\ < K (T) {I + \Zn\) Vs G [T„,r„+i] . (13) 

Since \xy^\ = |a;| \y\ for any x,y G (12) and (13) lead to 

l-Cn (s, cr)| < a: |cr| + a: (T) ^1 + Vs G [r„,T„+i], (14) 

where n = 0,..., N — 1 and £„ is as in (|^. Using Gronwall’s lemma, § and ( [l4| ) we deduce that 

(s)| < A:(T) ^1 + Vs G . (15) 

Decomposing 

an (t) an (t) - Pn (t) pI (f) 

as an (t) - Znzf^ - Zn (pn {t) - Zn)^ - {pn (t) - Z„) - (pn (t) - Zn) {pn (t) - Zn)^ WB have 

^n{t) = J Cn{s,an{s))ds - Zn(^j (H°/r„ (s) + 6° (s)) ds^ -j {B° Pn (s) + {s)) ds z^ 


j {B°n Un (s) + 6° (s)) ds(^J ( 5 ° Pn (s) + 6° (s)) ds^j , 


and so ( 13 ), ( 14 ) and ( 15 ) yields 


|ct„ (t)| < K (T) (^1 + (t-Tn) Vt G [r„,T„+i] 


( 16 ) 
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Iterating Q we obtain 

Zn+l ^ ^ (^n(s) f^n{s) («) + ^°(s) (s)) ds + Sn+1, 


where n = 0 ,..., N — 1 , nit) = max {n = 0 ,..., iV : < t} and 

n 

*^n+l ~ ^ ^ (^fc+l) M/c ("^fc+l) (’^fc+l) Vk ■ 


k—0 


Applying Holder’s inequality we get 

r'^rr + l 


/ (^°(.) M«(S) (s) + (s)) ds 

J to 

< (Tn+l - f M„(,) (s) + (s)) 

J in 


2g 


ds, 


and so (12) and ( 13 ) yield 


/ (Kis) l^n(s) (s) + (s)) ds 

J to 

Set Sq = 0 . For any n = 0 ,..., — 1 , 

e(| 5 „+i|") =E( 5 yi 5 „+i) 


2q 


<K{T)(l + 




'to 


i(s) I 


| 2 <? 


ds 


= X!® (^fc (cTfc (rfe+i) - (Tfc+i)^^ (rfc+i))?7fe) 

n d / r 

= (T-fc+l)^’^ + (m/c (Tfe+l)'^) 


fc=0 ^=1 


Since ak {jk+i) = E , ( 13 ) yields 

E 


(|5'n+l| ) ^ X/ (® I ) + ® iMfe ('Tfc+l)!^) < +00, 


fc=0 


(17) 


( 18 ) 


and so (S'n)„=o at is a (S^t„)„=o integrable martingale. According to the Burkholder-Davis- 

Gundy inequality we have 


E I max 
,k=0,. 


ax^ [si)"') < C,E = C,E ( ((ydfc in+i) 

’■■■’" 2 \fc=0 ^ 


rjk 


2 \ 1 


where Cq > 0 and stands for the j-th coordinate of the vector y. Applying Holder’s inequality yields 

(Tfc+l - (rfc+i - Tk)^^‘^ (j^y/dk {Tk+l)r]k'j ) / (Tfc+l - Tfc)^ 

/ " \ 9-1 n / _ 

- ^'^k+1 -Tk)j (i'*^+i “ I'fc) ( (\/crfc (Tfc+i) r?fc j j / (rfc+1 - Tkf 

\k=0 / k=0 ^ d 
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with 1/p + 1/g = 1. Using 


VO'fe ijk+l) 


Wk {Tk+i)\ we obtain 


E 


I max 



< {Tdy~^ 


< {Tdy~^ 




\/crk (Tfc+i) 


2q 


/c=0 


/c=0 


U, ^ E f (Tfe+l - Tfc) 


(rfc+1 - Tkf 

\^k ('T/c+l)! ^ 
'^k 


l%I^J 


Hence (16) yields 


E ^^max^ |5'fc|^‘^^ < K {T)E ^1 + ^ (rfc+i - r^) E( |zfc|^‘')^ . 

Using 0, @ and ( [I^ , together with Holder’s inequality, we get 

The discrete time Gronwall-Bellman lemma now leads to (101. 

We proceed to show 0. Using Holder’s inequality and (|l3|) we obtain 


(19) 


rr^+1 


(s) + 6° (s)) ds 


2q 


< {Tn+ 1 -T„f'‘ ^ J (|H°| |p„(s)| + |6° (s)|)^'^ds 

< K {T) {Tn+1 - Tnf'^ (l + \Zn\^'^'^ ■ 


By (16), 


\/ i'k'n+l) kjn 


2q 


< \J(Jn (r„+l) \rjn\‘^‘’ 

= l^n {'Tn+l)!'^ \Vuf‘^ 

< K (T) (^1 + (t„+i - TnY \r], 


1 1^ |29 


Hence 


E 


\/(In (Tn+l) T]n < K (T) (^1 + \Zn\^’^'^ (r„+i - Tn^ E ■ 


This implies (11), because 



1 (H° Hn (s) + &° (s)) ds 

+ 2^‘^-^E ( 

\J ("^n+l) 'Hn 

2q \ 

/?r„j . 


2q 




□ 
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Lemma 4 . 2 . Assume the hypothesis of Theorem \ 4 - 1 \ Let 

m 

Xn+I = f {tu, Zn) (t„+i “ t„) + ^ (t„, z„) . 


/c=l 


Then, for all n = 0 ,..., N — 1 , it is obtained that 


|E ((z„+l - - E(Xn+l/-i?T„)| < K (T) (t„+i - Tnf (1 + \Zn\) , 


(20) 


E 


and 


{{Zn+l - Zn) (Zn+l - 2n)^/?T„) “ E (Xn+ixl+l/5'r„) < K (T) (r„+i - T„)^ (^1 + , (21) 

( 22 ) 


E (^(2:„+i - Z„)^ (z„+i - Zn) {Zn+1 - Z„)^/5'r„) “ E (xl+lXn+lXl+l/-5r„) 
<iX(T) (T„+i-r„)" + 

Proof Since B° Zn + b°{T„) = f (r„, z„), 


(B° pin (s) + hn (s) - / ("^n, -^n)) ds 
{Bn (Lu (s) - ^^n) + K (s) - 6° (t„)) ds. 

Using Q and ( 13 ) we deduce that 

/ ■^TT+l 

{\pLn (s) - Z„| + S - T„) ds 

/ Tn + 1 pS 

J \B°pin (r) + 6° (s)| drds + K (r„+i - t„)^ 


Since 


< K (T) {Tn+l - Tn)'^ {I + \Zn\) ■ 


|E ((-2yi_|_i Zn) /^^Tn) E (Xn+l/^l5'r.i)| — \hn i'^n+l) Zn / ('^m ■^n) ("^n+l Ui)l 


( 23 ) 


@ yields 
From 

m 

E (Xn+lXn+l/-5r„) = / (t„, Zn) f (t„, Z„)^ (r„+i - T„)^ + X! 

we obtain 


/c=l 


IE (Xn+lXn+l/S^Tn) ^ ^ ('^n; ■ 2 ^n) ( 5 ^ (' 7 'mT^) 


fe=l 


< iX (T) (r„+i - T„)" (1 + |z„|') . ( 24 ) 


As in the proof of Lemma 4.1 we define (T„ (t) := an {t) — pin (t) pin (t) for any t G [t„, r„+i]. Then 


E 


(^{Zn+l - Zn) {Zn+1 “ •2n)^/5'r„) = {hn (Th+i) “ Zn) (hn {Tn+l) “ Zn)^ + an (t„+i) . 











Since 


cr„ (t„+i) = (T„ (t„+i) - ZnzJ - Z„ (^„ (t„+i) - Z„) - (^in (t„+i) - Z„) 

if-^n (^^n+l) -^n) (Mn (^n+l) -^n) 7 

applying (231 yields 

E (^{Zn+l - Zn) {Zn+1 - -Zn)^/5'r„) “ Cf-a (t„+i) + Z„zJ 

+ 2 :™/ ('Tn, (t„+i - t-^) + / (t„, z„) zJ (r„+i - t„) 

< if (T) (t„+i - T„)^ + |z„|^^ . 

Using (12), (13) and (14), together with Hypothesiswe deduce that 

|u„ (s, an (s)) - £n (r„, z„z^) | < if (T) (s - t„) (^1 + |z„|^^ , 


and so 


/ Tn+1 

\lLji (5, (7^1 (s)) ILfi {^n-) ) | 


< K (T) {Tn+l - TnY 


Therefore 


because 


m 

E ^(z„+i - Z„) (z„+i - Z„)^ -'^g^ (Tn, Zn) ("Tn, -^n))^ (Tn+l “ T„) 

/c=l 

< if (T) (r„+i - (^1 + |z„|^^ , 

m 

U„ (t„, z„z^) = z„/ (t„, z„)^ + / (t„, z„) z^ + X! ■ 

fe=i 


(25) 


Combining (24) with (25) we get (21). 

A careful computation shows 

E {xi+lXn+lXn+l/^r„) = f (r„, z„)^ / (r„, Zn) f (r„, z„)^ (t„+i - Tn)^ + f (r„, Zn)^ G„G^ (t„+i - T„)^ 

+/ (tti, z:n) (G„G^)^’ (r„+i - T„)^ + (G„G^) / (t„, z„)^ (r„+i - r„)^ , 

where G„ is the ’"-matrix whose (i,j)-th element is the i-th entry of (r„,z„). Similarly, 


E 


i,- 


(^{Zn+l - Znf (Z„+1 - Zn) (z„+i - Z„)^ /^Tr}j = ifJ-n (t„+i) - Z„)^ Cr„ (T-n-l-l) 

4” {g^n i'^n+l) Zn) an {^n+l) 

TCTyj (t 72 -|-i) (/Tn ('^n-t-l) Z^i) 

-I- (/i„ (t„+i) - z„)^ (^„ (Tn-i-i) - Z„) (/r„ (r„+i) - z„)^ . 
The last two inequalities imply (22), which completes the proof. □ 
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5 Numerical Simulations 


In this section, numerical simulations are presented in order to illustrate the performance of Scheme 
[2 This involves the numerical calculation of known expresions for functionals of two SDEs: a bilinear 
equation with random oscillatory dynamics, and a renowned nonlinear test equation. Fade method with 
scaling and squaring strategy (see, e.g., [T^]) was used to compute the exponential matrix in Q and (|^, 
whereas the squared root of the matrix (t„ (t„+i) — /i„ (r^i) ('Tn+i) in @ was computed by means 

of the singular value decomposition (see, e.g., i). vt in Q was set as a two-point distributed random 
variable with probability P{ri^ = ±1) = 1/2 for all n = 0,.., — 1 and fc = 1,.., m. All simulations were 

carried out in Matlab2014a. 


Example 1. Bilinear SDE with random oscillatory dynamics. 


dXt = a 


0 1 

-1 0 


Xtdt + Pi 


0 1 

-1 0 


XtdWl + P2 


1 0 
0 1 


XidWf 


(26) 


for all t G [0,12.5625], initial condition (Xq, Aq) = (1, 2), and parameters a = 10, pi = 0.1 and p 2 = 2pi. 
Since 


1 0 
0 1 


Xt = exp 


commutates with 


0 1 

-1 0 


{pi-pi)n 


, the solution of (|26|) is given by 

Wf + 


t + 


0 Pi 
-pi 0 


P2 0 
0 P2 


w: 


(27) 


(see, e.g., [T], p. 144). From Theorem 3 in [5], the mean mt and variance Vt of Xt are given by the 
expresions 

mt = Xq + L 2 exp(77t)uo 

and 

vec(yt) = Li exp(iJt)uo — vec{mtmj), 
where the matrices Li, L 2 , H and the vector uq are defined as 


(28) 

(29) 


H = 

■ A 

0 

0 

0 

0 ■ 
0 

e 

Mo = 

■ Mec(XoXT) ■ 
1 


0 

0 

C 



r 


Li — [ I/i O 4 ] € 


p4x8 


and 


T2 — [ O2 


x5 


O 2 


xl 


p2x8 


with 



1__ 

a 

2 

a 

2 

I—* to 

_1 


0 

a 

aXl 


■ 0 ■ 

A = 

—a 

P 2 

2 

-Pi 

2 

a 

e]R4x4^ C = 

—a 

0 

-aXl 

e and r = 

0 


—a 

-PI 

P 2 

a 


0 

0 

0 


1 


L Pi 

—a 

—a 

P2 \ 








First, we compare the exact values (28)-(29) for the mean and variance of Xt with their estimates 
obtained via Monte Carlo simulations. For this purpose, M realizations of the exact solution and 
Zn^ of the Schemej^were computed on an uniform time partition t„ = nA, with A = 1/2®, n = 0, ..,X, 
and N = 804. Then, with the estimates 

M 






d 


M 


and 


mr 


i=l 


M 

-is 


Ad 


i=l 
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for the mean, and 


M 


M 


i=l 


for the variance, the errors 



— nir 


and 

4 

411 = 

ml^ 


and 

41 = 

4!l = 

ml^ 


and 

4 = 

4!l = 



and 

4 = 

rol 

II 

14 ;^ 

-44 

and 

4 = 

4!! = 

14 ;^ 

-44 

and 

4 = 


1 


M 

M ^ 

i=l 




ml. — ml 


2 -^2 


vy-vy\ 


vy-vy\ 


vy-yy\ 


were evaluated. Here, for computing xll\ the realization of the Wiener process {Wl^,Wy was simulated 

n 

as ^ •\/AA/’(0, 1) for each k = 1,2, where A/’(0,1) is a Gaussian random 

f=i ^ ^ 

variable with zero mean and variance 1 . 

Figure 1 shows the exact values of versus their approximations mT„,Vr^ obtained from M = 

2^® simulations of Scheme Observe that there is not visual difference among these values. Table 
presents the errors = max{ei*j,} and = maxjerj,} of the estimated value of the mean and variance 


of (26 1 computed with different number of simulations M. As it was expected, these errors decrease as 
the number of simulations M increases. It is well known that the error e of the sampling mean of the 
Monte Carlo method decrease with the inverse of the square root of the number of simulations m, i-e., 

1 

e oc —— 

with 7 = 0.5. A roughly estimator of 7 for the errors and elrl was computed as minus the slope of 
the straight line fitted to the set of six points I (log 2 (Mfc), log 2 (eri (M/j))) : Mk = 2^, A: = 8,10,12,14,16, isj. 
Table shows the average 


N 


ZV] - 




for each type of error and its corresponding standard deviation 




\ 


1 

A^- 1 


N 




Results of Tables and together with Figure 1, indicate that the estimators for the mean and 
variance of (26 1 obtained by means of the simulations of the exact solution (27) and Scheme are quite 
similar. This is an expected result since the first two moments of the linear SDEs and Scheme 1 are 
’’equal” (up to the precision of the floating-point arithmetic in the numerical computation of the involved 
exponential and square root matrices). 
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/M 

2 ® 

210 

212 

214 

216 

to 

00 

elU 

0.10710 

0.05228 

0.04536 

0.01508 

0.00686 

0.00275 

et2l 

0.10643 

0.05025 

0.04469 

0.01433 

0.00647 

0.00304 

et3l 

0.43411 

0.25916 

0.18319 

0.29184 

0.07244 

0.03181 

et^l 

0.39102 

0.29529 

0.21413 

0.29496 

0.07726 

0.02753 

^5] 

0.23859 

0.14325 

0.15463 

0.16961 

0.05187 

0.02450 

eW 

0.27037 

0.02964 

0.02101 

0.02108 

0.01487 

0.00376 

e[2l 

0.27626 

0.04147 

0.02327 

0.02227 

0.01452 

0.00347 

e[3] 

0.92465 

0.35064 

0.18339 

0.15513 

0.06024 

0.02482 

e[41 

0.89503 

0.39518 

0.17646 

0.14678 

0.05655 

0.02346 

e[5] 

0.36642 

0.24892 

0.10664 

0.08899 

0.02785 

0.01101 


Table 1: Values of the errors and versus number of simulations M in the Example 1. 



eW 

el2l 

el3] 

eI4] 

ets] 

elU 

ePl 

e[3] 

e[4] 

e[5] 

7 

0.52 

0.53 

0.44 

0.44 

0.41 

0.44 

0.44 

0.44 

0.43 

0.45 

std 

0.16 

0.16 

0.20 

0.20 

0.21 

0.18 

0.19 

0.21 

0.21 

0.20 


Table 2: Average 7 and standard deviation std of the estimators for the rate of convergency 7 = 1/2 of the 
Monte Carlo simulations in the Example 1. 


/M 

2 ® 

210 

212 

214 

216 

21 ® 

rill 

0.0522 

0.0177 

0.0105 

0.0037 

0.0016 

0.0010 

r[2] 

0.0534 

0.0159 

0.0106 

0.0037 

0.0014 

0.0010 


Table 3: Relative error in the computation of the fnnctionals h\}^ and hr}, with different number of simulations 
M in the Example 1. 
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Figure 1: Integration of Example 1. Exact values of mt^,vt and their approximations fht^^Vt^ computed 
Monte Carlos with M = 2^® realizations of the SchemelH 


via 


In addition, let us compute the relative difference 

rW (M) = max {I )/h™ |} 

between the approximations 
1 


M 


li} = — > arctan 

i=l 


1 + 


Id 


and 


M 


/ijf] = — > arctan 

i=l 


1 + 




of the nonlinear functionals = E (arctan (l + (X^^)^)), with I = 1,2. Table displays the values of 
for different values of M. As it was also expected, rW goes to zero as the number of simulations M 
increases. Furthermore, Table shows that there is no significant difference between the estimates ob¬ 
tained from sampling the exact solution and Scheme]^ even though E (arctan (l -|- involves 

the computation of high order moments of Xr„ ■ 

The above simulation results illustrate the feasibility of Scheme for approximating functionals of 
linear SDEs with multiplicative noise. At this point is worth to mention that, with the uniform time 
partition consider here, the Euler scheme leads divergent results or computer overflows in the integration 


of the equation (26). 


Example 2. Nonautonomous nonlinear SDE 


r X} 1 


r -X? 1 


r 0 1 


cos(Xj+X^) 

. . 

— 

. XI _ 

dt -|- 

sin(Xl+Xf) 

VT+i 

dWf + 

Vi+t 

0 


dWf 


(30) 


with initial condition (Xq , Xq) = (1,1) and t G [0,10]. For this equation, E{(f){Xt)) = -|- log(l -I- 1), 

with (j){X) = \X'^. 
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Figure 2: Integration of Example 2. Exact Value: solid line. Schemejl] □ with A = 0.5, + with A = 0.25, 
* with A = 0.1. Euler with Romberg extrapolation: o with A = 0.05 and A = 0.1 


It is well-known from [T6j that via Monte Carlo simulations: 1) both, the Euler and the Milstein 
schemes with fixed stepsize A = 0.01 fail to approximate E{cj){Xt)); and 2) the second order method 
arising from Romberg’s extrapolation of the Euler scheme with stepsizes 0.02 and 0.01 gives a satisfactory 
approximation to E {(j) (Vj)), but fails when the stepsizes are 0.05 and 0.1. Similarly to the fourth figure 
in [T^, Figure 2 illustrates this last result for a Monte Carlo estimation with M = 10000 simulations. 

Figure 2 also shows the computation of E{(j){Xt)) via Monte Carlo method and Scheme]^ but on 
uniform time partitions with stepsizes A = 0.5,0.25,0.1 and M = 10000 simulations. In addition. Table 
[^provides the estimates e of the mean errors e = E{(f>{zj^)) — E{(j){Xx)) resulting from the integration 
of (30) via Scheme [^with different stepsizes. For this, the simulated trajectories i = 1, ...,K and 

j = 1,..., M, were are arranged into K = 100 batches of M = 10000 trajectories each for computing 


K 




with 


i=i 


1 

M 


M 

(At)). 


The 90% = 100(1 — a)% conhdence interval of the Student’s t distribution with K — 1 degrees for the 
mean error is given by 

[e — Ae, e -I- Ae] , 

where 

Ae = with ^1 = 

For comparison, the same estimate of the mean error for Euler scheme is also given in Table This 
illustrates again the better performance of the Scheme introduced in this paper. 
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<1 

1 

0.5 

0.25 

0.1 

Scheme |lj 

-2.2360 ±0.0093 

-0.4512 ±0.0067 

-0.0868 ± 0.0054 

0.0076 ± 0.0053 

Euler 

-2435.8 ± 1.7826 

-235.05 ±0.2192 

-32.031 ± 0.0361 

-5.7704 ±0.0101 


Table 4: Estimate e of the mean error E{(j)(zi\f)) — E{(f>{XT)) in the integration of (30l by means of Scheme]^ 
and the Euler scheme for different integration stepsizes A. 


6 Conclusions 

A weak Local Linearization scheme for stochastic differential equations with multiplicative noise was in¬ 
troduced. The scheme preserves the first two moments of the solution of linear SDEs and the mean square 
stability that such solution may have. The order-1 of weak convergence was proved and the practical 
performance of the scheme in the evaluation of functionals of linear and nonlinear SDEs was illustrated 
with numerical simulations. The simulations also showed the significant higher accuracy of the introduced 
scheme in comparison with the Euler scheme. 
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